      SUBROUTINE GALRKN(HW,M,ZZ,EGVL,FNORM,
     1              BEGVL,UMATL,VMATL,MX,MULT,A,B,H,HL,INT,INTH,
     2              MXSD2,FA,FB,FU,FV,FCS,FZ1,FZ2,FD,FE,
     3              DEPW,CGRADW,CINTW,RHOW,GRHOW,NPW,LHW,
     4              DEPB,CGRADB,CINTB,RHOB,GRHOB,NPB,LHB)
C
C     THIS SUBROUTINE FINDS COMPLEX EIGENVALUES EGVL AND EIGEN-
C     FUNCTIONS USING THE GALERKIN-MATRIX METHOD. THE GALERKIN 
C     MATRIX METHOD REDUCES THE PROBLEM TO A GENERALIZED MATRIX 
C     EIGENVALUE PROBLEM OF THE FORM: A UMATL = B UMATL DIAG(EGVL).
C     THE EIGENVALUES ARE THE SAME FOR BOTH PROBLEMS AND THE EIGEN-
C     FUNCTIONS MAY BE CONSTRUCTED FROM THE MATRIX EIGENVECTORS IN
C     UMATL. THE GENERALIZED MATRIX EIGENVALUE PROBLEM IS SYMMETRIZED 
C     USING THE CHOLESKY DECOMPOSITION OF B TO ASSURE OTHOGONALITY 
C     OF THE EIGENFUNCTIONS CONSTRUCTED FROM THE MATRIX EIGENVECTORS
C     IN UMATL.
C 
C     THIS VERSION ALLOWS PIECEWISE LINEAR VARIATONS IN THE SQUARE 
C     OF THE INDEX OF REFRACTION AND PIECEWISE EXPONENTIAL DENSITY    
C     VARIATION IN BOTH WATER AND BOTTOM. THE SQUARE OF THE INDEX 
C     OF REFRACTON CAN BE COMPLEX DUE TO ATTENUATION. 
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      COMMON /BLKEVN/ HB,CW,CB,FKW,FKB,ROHW,ROHB,ATEN
C
      DIMENSION INT(MX),INTH(MX)
      LOGICAL EIVEC
C
      REAL*8 ZZ(MX),ZERO,ESP,DREAL
C 
      DIMENSION DEPW(MX),RHOW(MX),GRHOW(MX)
      DIMENSION DEPB(MX),RHOB(MX),GRHOB(MX)
      REAL*8 B(MX,MX),VMATL(MX,MX)
      REAL*8 FA(MXSD2),FB(MXSD2),FU(MX,MX),FV(MX,MX)
      REAL*8 FCS(MX),FZ1(MX),FZ2(MX),FD(MX),FE(MX)
      COMPLEX*16 CGRADW(MX),CINTW(MX)
      COMPLEX*16 CGRADB(MX),CINTB(MX)
      COMPLEX*16 EGVL(MX),FNORM(MX),BEGVL(MX)
      COMPLEX*16 A(MX,MX),H(MX,MX),HL(MX,MX)
      COMPLEX*16 UMATL(MX,MX)
      COMPLEX*16 CDSRT,GAMA0
      COMPLEX*16 MULT(MX),CZERO,ONE
      COMPLEX*16 ALPHA1,BETA1
      COMPLEX*16 CKH,CSH
C
      ESP=.1D-30
      ZERO=DBLE(0.0)     
      ONE=DCMPLX(1.0,0.0)
      CZERO=DCMPLX(0.0,0.0)
C
C     FIND THE BASIC EIGENVALUES
C 
      CALL EGNVAL(HW,M,ZZ,MX)
C
C     FIND THE NORMALIZATION FACTORS FOR THE BASIC EIGENFUNCTIONS
C
      DO 100 I=1,M
C
      BEGVL(I)=DCMPLX(ZZ(I),ZERO)
      ALPHA1=CDSRT(FKW**2+FKB**2*(BEGVL(I)-ONE))
      BETA1=CDSRT(FKB**2*BEGVL(I))
      FNORM(I)=GAMA0(HW,ALPHA1,BETA1)
C
C      WRITE(6,50) I,BEGVL(I),FNORM(I),ALPHA1,BETA1
C      WRITE(9,50) I,BEGVL(I),FNORM(I),ALPHA1,BETA1
C   50 FORMAT(1X,I3,2X,4(2E12.4,1X))
C
  100 CONTINUE
C
C     COMPUTE THE REAL POSITIVE DEFINITE SYMMETRIC MATRIX B THAT
C     APPEARS ON THE RIGHT HAND SIDE OF THE GENERALIZED MATRIX EIGEN-
C     VALUE PROBLEM A UMATL = B UMATL DIAG(EGVL)
C
      CALL BMAT(HW,M,BEGVL,FNORM,B,MX,DEPW,RHOW,GRHOW,NPW,LHW,
     1          DEPB,RHOB,GRHOB,NPB,LHB)
C
C     PRINT THE MATRIX B
C
C      DO 170 I=1,M
C      DO 160 II=1,M
C      H(II,I)=DCMPLX(B(II,I),0.0D0)
C  160 CONTINUE
C  170 CONTINUE
C      CALL PRTMAT(M,M,'B   ',H,MX)
C
C     COMPUTE THE CHOLESKY DECOMPOSTION OF B
C
C                             T
C          B  =   VMATL  VMATL
C
C     INTO THE LOWER TRIANGULAR MATRIX VMATL AND ITS TRANSPOSE
C
      CALL CHOLDC(B,VMATL,M,MX)
C
C     CONSTRUCT THE COMPLEX SYMMETRIC MATRIX A THAT APPEARS ON THE
C     LEFT HAND SIDE OF THE GENERALIZED MATRIX EIGENVALUE PROBLEM
C     A UMATL = B UMATL DIAG(EGVL) 
C
      CALL AMAT(HW,M,BEGVL,FNORM,A,MX,
     1          DEPW,CGRADW,CINTW,RHOW,GRHOW,NPW,LHW,
     2          DEPB,CGRADB,CINTB,RHOB,GRHOB,NPB,LHB)
C
C     PRINT THE MATRIX A
C
C      CALL PRTMAT(M,M,'A   ',A,MX)
C
C     USE FORWARD SUBSTITUTION TWICE TO COMPUTE THE COMPLEX SYMMETRIC 
C     MATRIX
C                -1         T -1
C           VMATL   A  VMATL     ,    
C
C     WHICH IS STORED IN A, AND WHOSE EIGENVALUES ARE THE SAME AS THE
C     GENERALIZED EIGENVALUE PROBLEM A UMATL = B UMATL DIAG(EGVL) 
C
      CALL FORESUB(VMATL,HL,A,M,MX)
C
C     DEFINE THE TRANSPOSE OF HL
C
      DO 230 I=1,M
      DO 210 II=1,M
      H(II,I)=HL(I,II)
  210 CONTINUE
  230 CONTINUE
C
C     SOLVE FOR THE SYMMETRIIC MATIX
C
      CALL FORESUB(VMATL,A,H,M,MX)
C
C     PRINT THE SYMMETRIC MATRIX A 
C
C      CALL PRTMAT(M,M,'SA  ',A,MX)
C
C     STORE THE REAL AND IMAGINARY PARTS OF THE DIAGONAL AND SUPER
C     DIAGONAL OF THE SYMMETRIC MATRIX A, BY COLUNMS, TO BE PASSED 
C     TO THE JACOBI TYPE ALGORITHM CSYMEIG
C
      DO 270 J=1,M
      DO 250 I=1,J
      JJI=(J*(J-1)/2)+I
      FA(JJI)=DREAL(A(I,J))
      FB(JJI)=DIMAG(A(I,J))
  250 CONTINUE
  270 CONTINUE
      EIVEC=.TRUE.
C
C     FIND THE EIGENVALUES AND EIGENVECTORS USING THE JACOBI TYPE 
C     ALGORITHM. 
C
      CALL CSYMEIG(FA,FB,FU,FV,FCS,FZ1,FZ2,FD,FE,M,EIVEC,MX,MXSD2)
C
C     THE REAL AND IMAGINARY PARTS OF THE EIGENVALUES, FROM THE
C     DIAGONAL ENTRIES OF FA AND FB, ARE STORED IN FD AND FD.
C
      DO 300 I=1,M
      II=I*(I+1)/2
      FD(I)=FA(II)
      FE(I)=FB(II)
      INT(I)=I
  300 CONTINUE
C
C     INDEX THE EIGENVALUES BY INCREASING REAL PART
C
      CALL INDEXX(M,FD,INT,INTH,MX)
C
C     DEFINE THE COMPLEX EIGENVALUES IN ORDER BY INCREASING
C     REAL PART
C
C      WRITE(9,*) ' '
C      WRITE(9,*) 'HORIZONTAL WAVE NUMBERS:'
C      WRITE(9,*) ' '
C
      DO 400 I=1,M
c
c      if(i .ne. int(M-i+1)) then
c      write(6,*) i,int(M-i+1)
c      write(*,*) i,int(M-i+1)
c      end if
c
      EGVL(I)=DCMPLX(FD(INT(I)),FE(INT(I)))
C
C      CKH=FKB*CDSRT(ONE-EGVL(I))
C      WRITE(6,390) I,EGVL(I)
C      WRITE(9,390) I,CKH
C  390 FORMAT(1X,I5,2G23.14)
C
  400 CONTINUE
C
C      WRITE(9,*) ' '
C      WRITE(9,*) 'HORIZONTAL SLOWNESS:'
C      WRITE(9,*) ' '
C      OMEGA=CB*FKB
C      DO 410 I=1,M
C      CKH=FKB*CDSRT(ONE-EGVL(I))
C      CSH=1000.0D0*CKH/OMEGA
C      WRITE(6,390) I,CSH
C      WRITE(9,390) I,CSH
C  410 CONTINUE
C
C     DEFINE THE COMPLEX EIGENVECTORS
C
      DO 500 I=1,M
      DO 490 II=1,M
      UMATL(II,I)=DCMPLX(FU(II,INT(I)),FV(II,INT(I)))
  490 CONTINUE
  500 CONTINUE
C
C     PRINT THE MATRIX OF EIGENVECTORS
C
C      CALL PRTMAT(M,M,'UMAT',UMATL,MX)
C
C     IN SOME CASES, WHEN THERE IS LOW ATTENUATION, NON-PHYSICAL
C     EIGENVALUES CAN BE F0UND. IN COUPLE, PHYSICAL EIGNEVALUES
C     HAVE NEGATIVE IMAGINARY PARTS (THIS IS BEFORE THEY ARE
C     SUBTRACTED FROM ONE AND A SQUARE ROOT, WITH A POSITIVE
C     IMAGINARY PART, IS TAKEN TO GET WAVE NUMBERS).
C
      DO 600 I=1,M
C
      IF(DIMAG(EGVL(I)) .GT. ESP) THEN
C
      WRITE(9,*) '*** NON-PHYSICAL EIGENVLAUE FOUND FOR',
     1           ' MODE= ',I,' IN GALRKN ***'
      WRITE(9,*) 'Try putting a small amount of attenuation ',
     1           'in the water layer (1e-6 dB/wl)'
      WRITE(6,*) '*** NON-PHYSICAL EIGENVLAUE FOUND FOR',
     1           ' MODE= ',I,' IN GALRKN ***'
      WRITE(6,*) 'Try putting a small amount of attenuation ',
     1           'in the water layer (1e-6 dB/wl)'
C
      IFLAG=1
      CALL ABORTC(IFLAG)
      END IF
C
  600 CONTINUE
C
C     DEFINE THE TRANSPOSE OF VMATL
C
      DO 800 I=1,M
      DO 700 II=1,M
      H(II,I)=DCMPLX(VMATL(I,II),0.0D0)
  700 CONTINUE
  800 CONTINUE
C
C     RECOVER THE EIGENVECTORS OF THE GENERALIZED MATRIX EIGENVALUE
C     PROBLEM
C                 T -1  
C            VMATL     UMATL    
C
C     USING BACKSUBSTITUTION AND STORE THE RESULTS IN HL.
C
      CALL BACKSUB(H,HL,UMATL,M,MX)
C
C     RETURN THE EIGENVECTORS IN UMATL AND STORE THE REAL SYMMETRIC
C     POSITIVE DEFINITE MATRIX B IN THE COMPLEX MATRIX A FOR ORTHO-
C     GONALITY CHECK
C
      DO 1000 I=1,M
      DO 900 II=1,M
      UMATL(II,I)=HL(II,I)
      A(II,I)=DCMPLX(B(II,I),0.0D0)
  900 CONTINUE
 1000 CONTINUE
C
C     CHECK ORTHOGONALITY
C
C      CALL TMMULT(M,UMATL,A,UMATL,MX,H)
C      CALL PRTMAT(M,M,'ORTH',A,MX)
C
      RETURN 
      END      
